clear all

use "$root\data\use.dta", replace


foreach crop in corn soybeans {
	local c_name = strproper("`crop'")

	preserve
	local c = cond("`crop'" == "corn", "c", "s")
	rename `c'_d*_* d*_*
	rename `c'_tot_* tot_*
	rename `c'_pdsi_w_* pdsi_w_*
	rename `c'_pmdi_w_* pmdi_w_*
	rename `c'_tbin*_* tbin*_*
	rename `c'_prec_* prec_*
	rename `crop'_production00 production00

**## Planned planting acres
	reghdfe ln`crop'_pplanted d*_b d*_p if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100 , absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_pplanted_wo
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p

	reghdfe ln`crop'_pplanted d*_b d*_p tbin*_b prec_b* tbin*_p prec_p* if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_pplanted_w
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p
	
	reghdfe ln`crop'_pplanted d*_b d*_p pdsi_w_b tbin*b pdsi_w_p tbin*p if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_pplanted_wet
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p

	
**## Planting ratio
	reghdfe `crop'_plratio d*_b d*_p if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100 , absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_plratio_wo
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p

	reghdfe `crop'_plratio d*_b d*_p tbin*_b prec_b* tbin*_p prec_p* if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_plratio_w
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p
	
	reghdfe `crop'_plratio d*_b d*_p pdsi_w_b tbin*b pdsi_w_p tbin*p if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_plratio_wet
	test d0_b d1_b d2_b d3_b d0_p d1_p d2_p d3_p

	
**## Harvest ratio
	reghdfe `crop'_hratio d*_b d*_p d*_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_wo

	reghdfe `crop'_hratio d*_b d*_p d*_g tbin*_b prec*_b* tbin*_p prec*_p* tbin*_g prec*_g* if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_w
	
	reghdfe `crop'_hratio d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_wet
	
	reghdfe `crop'_hratio d*_g pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_hratio_omit


**## Yield
	reghdfe ln`crop'_yield d*_b d*_p d*_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_wo

	reghdfe ln`crop'_yield d*_b d*_p d*_g tbin*_b prec*_b* tbin*_p prec*_p* tbin*_g prec*_g* if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_w

	reghdfe ln`crop'_yield d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_yield_wet
	
		
**## Production
	reghdfe ln`crop'_production d*_b d*_p d*_g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_wo

	reghdfe ln`crop'_production d*_b d*_p d*_g tbin*_b prec*_b* tbin*_p prec*_p* tbin*_g prec*_g* if  irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_w

	reghdfe ln`crop'_production d*_b d*_p d*_g pdsi_w_b tbin*b pdsi_w_p tbin*p pdsi_w_g tbin*g if irrigated_ratio<0.1 & N_year_`crop'>10 & lng>-100, absorb(fips i.state#c.year i.state#c.year2 i.year) vce(cluster fips stateyear)
	estimate store `crop'_production_wet

restore
}


**## Plots
coefplot (corn_pplanted_wo, label(No Weather Controls) offset(-0.15) ciopts(color(black) )) /// 
	(corn_pplanted_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) )  msymbol(T)) ///
	(corn_pplanted_wet, label(Temp. and Wetness Controls) offset(0.15) ciopts(color(black)) msymbol(S)), bylabel(Corn) || ///
	(soybeans_pplanted_wo, label(No Weather Controls) offset(-0.15) ciopts(color(black)  )) /// 
	(soybeans_pplanted_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) ) msymbol(T)) /// 
	(soybeans_pplanted_wet, label(Temp. and Wetness Controls) offset(0.15) ciopts(color(black)) msymbol(S)), bylabel(Soybeans) ||, ///
	mcolor(black) mfcolor(white) vertical keep(d*_b d*_p) yline(0, lcolor(gray) lwidth(0.2)) /// 
	byopts(compact cols(1) graphregion(fcolor(white))) subtitle(, fcolor(white) bmargin(top) bcolor(white)) ///
	groups(d*_b = Pre-planting d*p=Planting) ytitle("ln(Planned Acres)") ylab(,labs(small)) xlabel(1 "D0" 2 "D1" 3 "D2" 4 "D3" 6 "D0" 7 "D1" 8 "D2" 9 "D3")
graph export "$root/results_replication/figures/f1_main_pplanted.png", as(png) name("Graph")  width(3200) height(2400) replace

coefplot (corn_plratio_wo, label(No Weather Controls) offset(-0.15) ciopts(color(black) )) /// 
	(corn_plratio_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) )  msymbol(T)) ///
	(corn_plratio_wet, label(Temp. and Wetness Controls) offset(0.15) ciopts(color(black)) msymbol(S)), bylabel(Corn) || ///
	(soybeans_plratio_wo, label(No Weather Controls) offset(-0.15) ciopts(color(black)  )) /// 
	(soybeans_plratio_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) ) msymbol(T)) /// 
	(soybeans_plratio_wet, label(Temp. and Wetness Controls) offset(0.15) ciopts(color(black)) msymbol(S)), bylabel(Soybeans) ||, ///
	mcolor(black) mfcolor(white) vertical keep(d*_b d*_p) yline(0, lcolor(gray) lwidth(0.2)) /// 
	byopts(compact cols(1) graphregion(fcolor(white))) subtitle(, fcolor(white) bmargin(top) bcolor(white)) ///
	groups(d*_b = Pre-planting d*_p = Planting) ytitle("Planted Ratio") ylab(,labs(small)) xlabel(1 "D0" 2 "D1" 3 "D2" 4 "D3" 6 "D0" 7 "D1" 8 "D2" 9 "D3")
graph export "$root/results_replication/figures/f2_main_plratio.png", as(png) name("Graph")  width(3200) height(2400) replace

coefplot (corn_hratio_wo, label(No Weather Controls) offset(-0.22) ciopts(color(black) )) /// 
	(corn_hratio_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) )  msymbol(T)) ///
	(corn_hratio_wet, label(Temp. and Wetness Controls) offset(0.22) ciopts(color(black)) msymbol(S)), bylabel(Corn) || ///
	(soybeans_hratio_wo, label(No Weather Controls) offset(-0.22) ciopts(color(black)  )) /// 
	(soybeans_hratio_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) ) msymbol(T)) /// 
	(soybeans_hratio_wet, label(Temp. and Wetness Controls) offset(0.22) ciopts(color(black)) msymbol(S)), bylabel(Soybeans) ||, ///
	mcolor(black) mfcolor(white) vertical keep(d*_b d*_p d*_g) yline(0, lcolor(gray) lwidth(0.2)) /// 
	byopts(compact cols(1) graphregion(fcolor(white))) subtitle(, fcolor(white) bmargin(top) bcolor(white)) ///
	groups(d*_b = Pre-planting d*_p = Planting d*_g = Growing) ytitle("Harvested Ratio") ylab(,labs(small)) xlabel(1 "D0" 2 "D1" 3 "D2" 4 "D3" 6 "D0" 7 "D1" 8 "D2" 9 "D3" 11 "D0" 12 "D1" 13 "D2" 14 "D3")
graph export "$root/results_replication/figures/f3_main_hratio.png", as(png) name("Graph")  width(3200) height(2400) replace

foreach outcome in yield production {
    local y_name = strproper("`outcome'")
	local order = cond("`outcome'" == "hratio", 3, cond("`outcome'" == "yield", 4, 5))
	coefplot (corn_`outcome'_wo, label(No Weather Controls) offset(-0.22) ciopts(color(black) )) /// 
	(corn_`outcome'_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) )  msymbol(T)) ///
	(corn_`outcome'_wet, label(Temp. and Wetness Controls) offset(0.22) ciopts(color(black)) msymbol(S)), bylabel(Corn) || ///
	(soybeans_`outcome'_wo, label(No Weather Controls) offset(-0.22) ciopts(color(black)  )) /// 
	(soybeans_`outcome'_w, label(Temp. and Precip. Controls) offset(0) ciopts(color(black) ) msymbol(T)) /// 
	(soybeans_`outcome'_wet, label(Temp. and Wetness Controls) offset(0.22) ciopts(color(black)) msymbol(S)), bylabel(Soybeans) ||, ///
	mcolor(black) mfcolor(white) vertical keep(d*_b d*_p d*_g) yline(0, lcolor(gray) lwidth(0.2)) /// 
	byopts(compact cols(1) graphregion(fcolor(white))) subtitle(, fcolor(white) bmargin(top) bcolor(white)) ///
	groups(d*_b = Pre-planting d*_p = Planting d*_g = Growing) ytitle("ln(`y_name')") ylab(,labs(small)) xlabel(1 "D0" 2 "D1" 3 "D2" 4 "D3" 6 "D0" 7 "D1" 8 "D2" 9 "D3" 11 "D0" 12 "D1" 13 "D2" 14 "D3")
	graph export "$root/results_replication/figures/f`order'_main_`outcome'.png", as(png) name("Graph")  width(3200) height(2400) replace
}